
# delimit;
clear;
clear matrix;
capture log close;
capture matrix close;
set more off;
set mem 400m;



cd D:\Data\WorkData\703201\vibeke\GP;
log using "Happy Doctor_graphs_log.txt", text replace;


*****************************************************;
* FIGURE 1*******************************************;
*****************************************************;

* log birth weight;

use infant2.dta,clear;
drop if zone==1;
drop if pre1==.;
keep if kpnrm==1 | kpnrm==2;
keep if oob==1;
char _dta[omit] "prevalent";
set seed 210177;
keep if agem<27;

for Z in any 1 2:
\xi: xtreg lgbw yb4 yb6-yb9 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  male amx1 amx3 amx4 singl west nowest udlm qm2-qm4 o4m udlf qf2-qf4 mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if kpnrm==Z, i(kompnrm)  fe robust 
\predict xZ if e(sample);
gen bb=x1;
replace bb=x2 if kpnrm==2;
drop x1 x2;
gen double half= hofd(fodtdato);
format half %th;
egen mbb=mean(bb), by (kpnrm half );
sort half kpnrm;
keep if kpnrm!=kpnrm[_n-1];
line mbb half if kpnrm==1, xaxis(1 2) yaxis(1) c(l) clcolor(black) clpat(solid) clwidth(medthick)
||line mbb half if kpnrm==2 , xaxis (1 2) yaxis(1) c(l) clcolor(black) clpat(longdash) clwidth(medthick)
xline(55.5 57)
xlabel("",axis(2)) xtitle("",axis(2))
ylabel(8.00 (0.02) 8.15, angle(horizontal))

ytitle("Log points", axis (1))
xtitle("Birth halfyear")
legend(label(1 "Treatment group") label(2 "Control group"))
title("Log birth weight, Mother's age<27")
;
graph save FIG1_LGBW_YOUNG,replace;
********************************';


*>30;
use infant2.dta,clear;
drop if zone==1;
drop if pre1==.;
keep if kpnrm==1 | kpnrm==2;
keep if oob==1;
char _dta[omit] "prevalent";
set seed 210177;
keep if agem>30;
for Z in any 1 2:
\xi: xtreg lgbw yb4 yb6-yb9 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  male amx1 amx3 amx4 singl west nowest udlm qm2-qm4 o4m udlf qf2-qf4 mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if kpnrm==Z, i(kompnrm)  fe robust 
\predict xZ if e(sample);
gen bb=x1;
replace bb=x2 if kpnrm==2;
drop x1 x2;
gen double half= hofd(fodtdato);
format half %th;
egen mbb=mean(bb), by (kpnrm half );
sort half kpnrm;
keep if kpnrm!=kpnrm[_n-1];
line mbb half if kpnrm==1, xaxis(1 2) yaxis(1) c(l) clcolor(black) clpat(solid) clwidth(medthick)
||line mbb half if kpnrm==2 , xaxis (1 2) yaxis(1) c(l) clcolor(black) clpat(longdash) clwidth(medthick)
xline(55.5 57)
xlabel("",axis(2)) xtitle("",axis(2))
ytitle("Log points", axis (1))
ylabel(8 (0.02) 8.15, angle(horizontal))
xtitle("Birth halfyear")
legend(label(1 "Treatment group") label(2 "Control group"))
title("Log birth weight, Mother's age>30")
;
graph save FIG1_LGBW_OLD,replace;
******************************************************;
* preterm birth;
*<27;

use infant2.dta,clear;
drop if pre1==.;
keep if yob>1984;

drop if zone==1;
keep if kpnrm==1 | kpnrm==2;
keep if oob==1;
char _dta[omit] "prevalent";
set seed 210177;
keep if agem<27;
generate float nydat = qofd(fodtdato);

for Z in any 1 2:
\xi: xtreg ptb yb4 yb6-yb9 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  male amx1 amx3 amx4 singl west nowest udlm qm2-qm4 o4m udlf qf2-qf4 mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if kpnrm==Z, i(kompnrm)  fe robust 
\predict xZ if e(sample);
gen bb=x1;
replace bb=x2 if kpnrm==2;
drop x1 x2;
gen double half= hofd(fodtdato);
format half %th;
egen mbb=mean(bb), by (kpnrm half );
sort half kpnrm;
keep if kpnrm!=kpnrm[_n-1];
line mbb half if kpnrm==1, xaxis(1 2) yaxis(1) c(l) clcolor(black) clpat(solid) clwidth(medthick)
||line mbb half if kpnrm==2 , xaxis (1 2) yaxis(1) c(l) clcolor(black) clpat(longdash) clwidth(medthick)
xline(55.5 57)
xlabel("",axis(2)) xtitle("",axis(2))
ytitle("frequency", axis (1))
ylabel(0.04 (0.02) 0.16, angle(horizontal))
xtitle("Birth halfyear")
legend(label(1 "Treatment group") label(2 "Control group"))
title("Preterm birth, Mother's age<27")
;
graph save FIG1_PTB_YOUNG,replace;

********************************';
*>30;
use infant2.dta,clear;
gen llbw= (bw<2000);
drop if pre1==.;
drop if zone==1;
keep if kpnrm==1 | kpnrm==2;
keep if oob==1;
char _dta[omit] "prevalent";
set seed 210177;
keep if agem>30;

generate float nydat = qofd(fodtdato);

for Z in any 1 2:
\xi: xtreg ptb yb4 yb6-yb9 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  male amx1 amx3 amx4 singl west nowest udlm qm2-qm4 o4m udlf qf2-qf4 mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if kpnrm==Z, i(kompnrm)  fe robust 
\predict xZ if e(sample);
gen bb=x1;
replace bb=x2 if kpnrm==2;
drop x1 x2;
gen double half= hofd(fodtdato);
format half %th;
egen mbb=mean(bb), by (kpnrm half );
sort half kpnrm;
keep if kpnrm!=kpnrm[_n-1];
line mbb half if kpnrm==1, xaxis(1 2) yaxis(1) c(l) clcolor(black) clpat(solid) clwidth(medthick)
||line mbb half if kpnrm==2 , xaxis (1 2) yaxis(1) c(l) clcolor(black) clpat(longdash) clwidth(medthick)
xline(55.5 57)
xlabel("",axis(2)) xtitle("",axis(2))
ytitle("frequency", axis (1))
ylabel(0.04 (0.02) 0.16, angle(horizontal))
xtitle("Birth halfyear")
legend(label(1 "Treatment group") label(2 "Control group"))
title("Preterm birth, Mother's age>30")
;
graph save FIG1_LGBW_OLD,replace;

*********************************************;
* very preterm birth;
*<27;

use infant2.dta,clear;
drop if pre1==.;
keep if yob>1984;
gen llbw= (bw<2000);
drop if zone==1;
keep if kpnrm==1 | kpnrm==2;
keep if oob==1;
char _dta[omit] "prevalent";
set seed 210177;
keep if agem<27;
generate float nydat = qofd(fodtdato);

for Z in any 1 2:
\xi: xtreg vptb yb4 yb6-yb9 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  male amx1 amx3 amx4 singl west nowest udlm qm2-qm4 o4m udlf qf2-qf4 mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if kpnrm==Z, i(kompnrm)  fe robust 
\predict xZ if e(sample);
gen bb=x1;
replace bb=x2 if kpnrm==2;
drop x1 x2;
gen double half= hofd(fodtdato);
format half %th;
egen mbb=mean(bb), by (kpnrm half );
sort half kpnrm;
keep if kpnrm!=kpnrm[_n-1];
line mbb half if kpnrm==1, xaxis(1 2) yaxis(1) c(l) clcolor(black) clpat(solid) clwidth(medthick)
||line mbb half if kpnrm==2 , xaxis (1 2) yaxis(1) c(l) clcolor(black) clpat(longdash) clwidth(medthick)
xline(55.5 57)
xlabel("",axis(2)) xtitle("",axis(2))
ytitle("frequency", axis (1))
ylabel(0.00 (0.01) 0.06, angle(horizontal))
xtitle("Birth halfyear")
legend(label(1 "Treatment group") label(2 "Control group"))
title("Very preterm birth, Mother's age<27")
;
graph save FIG1_VPTB_YOUNG,replace;

********************************';

*>30;
use infant2.dta,clear;
drop if pre1==.;
drop if zone==1;
keep if kpnrm==1 | kpnrm==2;
keep if oob==1;
char _dta[omit] "prevalent";
set seed 210177;
keep if agem>30;
generate float nydat = qofd(fodtdato);


for Z in any 1 2:
\xi: xtreg vptb yb4 yb6-yb9 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  male amx1 amx3 amx4 singl west nowest udlm qm2-qm4 o4m udlf qf2-qf4 mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if kpnrm==Z, i(kompnrm)  fe robust 
\predict xZ if e(sample);
gen bb=x1;
replace bb=x2 if kpnrm==2;
drop x1 x2;
gen double half= hofd(fodtdato);
format half %th;
egen mbb=mean(bb), by (kpnrm half );
sort half kpnrm;
keep if kpnrm!=kpnrm[_n-1];
line mbb half if kpnrm==1, xaxis(1 2) yaxis(1) c(l) clcolor(black) clpat(solid) clwidth(medthick)
||line mbb half if kpnrm==2 , xaxis (1 2) yaxis(1) c(l) clcolor(black) clpat(longdash) clwidth(medthick)
xline(55.5 57)
xlabel("",axis(2)) xtitle("",axis(2))
ytitle("frequency", axis (1))
xtitle("Birth halfyear")
ylabel(0.00 (0.01) 0.06, angle(horizontal))
legend(label(1 "Treatment group") label(2 "Control group"))
title("Very preterm birth, Mother's age>30")
;
graph save FIG1_VPTB_OLD,replace;


* fetal growth;
*<27;

use infant2.dta,clear;
drop if pre1==.;
keep if yob>1984;
gen llbw= (bw<2000);
drop if zone==1;
keep if kpnrm==1 | kpnrm==2;
keep if oob==1;
char _dta[omit] "prevalent";
set seed 210177;
keep if agem<27;
generate float nydat = qofd(fodtdato);


for Z in any 1 2:
\xi: xtreg fg yb4 yb6-yb9 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  male amx1 amx3 amx4 singl west nowest udlm qm2-qm4 o4m udlf qf2-qf4 mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if kpnrm==Z, i(kompnrm)  fe robust 
\predict xZ if e(sample);
gen bb=x1;
replace bb=x2 if kpnrm==2;
drop x1 x2;
gen double half= hofd(fodtdato);
format half %th;
egen mbb=mean(bb), by (kpnrm half );
sort half kpnrm;
keep if kpnrm!=kpnrm[_n-1];
line mbb half if kpnrm==1, xaxis(1 2) yaxis(1) c(l) clcolor(black) clpat(solid) clwidth(medthick)
||line mbb half if kpnrm==2 , xaxis (1 2) yaxis(1) c(l) clcolor(black) clpat(longdash) clwidth(medthick)
xline(55.5 57)
xlabel("",axis(2)) xtitle("",axis(2))
ytitle("grams per week", axis (1))
xtitle("Birth halfyear")
ylabel(80 (2) 90, angle(horizontal))
legend(label(1 "Treatment group") label(2 "Control group"))
title("Fetal growth, Mother's age<27")
;
graph save FIG1_FG_YOUNG.gph,replace;

********************************';
*>30;
use infant2.dta,clear;
drop if pre1==.;
drop if zone==1;
keep if kpnrm==1 | kpnrm==2;
keep if oob==1;
char _dta[omit] "prevalent";
set seed 210177;
keep if agem>30;
generate float nydat = qofd(fodtdato);


for Z in any 1 2:
\xi: xtreg fg yb4 yb6-yb9 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  male amx1 amx3 amx4 singl west nowest udlm qm2-qm4 o4m udlf qf2-qf4 mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if kpnrm==Z, i(kompnrm)  fe robust 
\predict xZ if e(sample);
gen bb=x1;
replace bb=x2 if kpnrm==2;
drop x1 x2;
gen double half= hofd(fodtdato);
format half %th;
egen mbb=mean(bb), by (kpnrm half );
sort half kpnrm;
keep if kpnrm!=kpnrm[_n-1];
line mbb half if kpnrm==1, xaxis(1 2) yaxis(1) c(l) clcolor(black) clpat(solid) clwidth(medthick)
||line mbb half if kpnrm==2 , xaxis (1 2) yaxis(1) c(l) clcolor(black) clpat(longdash) clwidth(medthick)
xline(55.5 57)
xlabel("",axis(2)) xtitle("",axis(2))
ytitle("grams per week", axis (1))
xtitle("Birth halfyear")
ylabel(80 (2) 90, angle(horizontal))
legend(label(1 "Treatment group") label(2 "Control group"))
title("Fetal growth, Mother's age>30")
;
graph save FIG1_FG_OLD,replace;


**********************************************************;

*************************************************************************************************************************;
** APPENDIX FIGURE A1: REGIONAL TAXATION IN THE TREATMENT AND  CONTROL AREA;
*** DATA FROM STATISTIC DENMARKS DYNAMIC DATABASES: STATBANK AT WWW.DST.DK 




*************************************************************************************************************************;
*RECODNING FOR APPENDIX FIGURE A2 & A3 FIGURES BY AGE GROUPS;
*************************************************************************************************************************;
use infant2.dta,clear;
drop if zone==1;
drop if pre1==.;
keep if kpnrm==1 | kpnrm==2;
*keep if oob==1;
char _dta[omit] "prevalent";
set seed 210177;

gen age=.;
replace age=20 if agem<=20;
replace age=22 if agem>20 & agem<=22;
replace age=24 if agem>22 & agem<=24;
replace age=26 if agem>24 & agem<=26;
replace age=28 if agem>26 & agem<=28;
replace age=30 if agem>28 & agem<=30;
replace age=32 if agem>30 & agem<=32;
replace age=34 if agem>32 & agem<=34;
replace age=36 if agem>34 & agem<=36;
replace age=38 if agem>36 & agem<=38;
replace age=40 if agem>38 ;

label define age  20 "-20" 22 "21-22" 24 "23-24" 26 "25-26" 28 "27-28" 30 "29-30" 32 "31-32" 34 "33-34" 36 "35-36" 38 "37-38" 40 "39+";
label values age age;
label variable age "age";
save tmp.dta,replace;




****************************************************************************************************************;
*APPENDIX FIGURE A2: THE DISTRIBUTION OF LOW BIRTH WEIGHT AND PRETERM BIRTH FOR FIST TIME MOTHERS, BY AGE ;
****************************************************************************************************************;


graph bar (mean) lbw ptb if oob==1,  legend( label(1 "Low birth weight") label(2 "Preterm birth") ) over(age) ytitle("Probability") ;
graph save FIGA2.gph,replace;
# delimit ;


****************************************************************************************************************;
* APPENDIX FIGURE A3.--THE EFFECT OF CAPITATION CONTRACTS ON LOG BIRTH WEIGHT AND PRETERM BIRTH, BY MATERNAL AGE;
****************************************************************************************************************;

use tmp.dta,clear;
capture program drop  keepmat;
prog define keepmat;
matrix C=e(b);
matrix V=e(V);
scalar W=sqrt(V[1,1]);
scalar y=(C[1,1]/W);
scalar O=sqrt(V[2,2]);
scalar p=(C[1,2]/O);
matrix F=F\(C[1,1],W,y, e(N), C[1,2],O,p);
end;
matrix F=(0,0,0,0,0,0,0);

for A in num 16/40:

\quietly xi: xtreg lgbw preca pre case  yb1-yb4 yb6 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  o4m  male  amx2 amx3 amx4 singl west nowest udlm qm2-qm4 udlf mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if agem>=A & agem<=(A+6), i(kompnrm)  fe robust 
\keepmat
;



matrix list F;




use tmp.dta,clear;
capture program drop  keepmat;
prog define keepmat;
matrix C=e(b);
matrix V=e(V);
scalar W=sqrt(V[1,1]);
scalar y=(C[1,1]/W);
scalar O=sqrt(V[2,2]);
scalar p=(C[1,2]/O);
matrix G=G\(C[1,1],W,y, e(N), C[1,2],O,p);
end;
matrix G=(0,0,0,0,0,0,0);
for A in num 16/40:
\quietly: xi: xtreg ptb preca pre case  yb1-yb4 yb6 mb1-mb8 mb10-mb12 Docsex Docage Docagesq  o4m  male  amx2 amx3 amx4 singl west nowest udlm qm2-qm4 udlf mpnrf mo4f mo4m  msingl mudlf mudlm mqincm mqincf if agem>=A & agem<=(A+6), i(kompnrm)  fe robust 
\keepmat
;
matrix list G;

*************************************************************************************************;
* obs coding by hand;
*1: save the two matrix lists F & G as stata-files (I did a cut and paste); 
*2: keep first column (the mean);
*3 recode row number to age according to the loop 16/42;
* 4 save the data as the dataset age.dta;
*************************************************************************************************;




use age.dta,clear;
destring, replace;
*rename (var1 var2 var3) (age lgbw ptb);
line lgbw age,
ylabel(-0.08(0.01)0.01, axis(1) gmin angle(horizontal)) 
 xlabel(16(2)40, axis(1))
ytitle("log points")
xtitle("maternal age + 6 years")
subtitle("THE EFFECT ON LOG BIRTH WEIGHT");
graph save FIGA3_lbw.gph,replace;

line ptb age,
ylabel(0(0.01)0.08, axis(1) gmin angle(horizontal)) 
 xlabel(16(2)40, axis(1))
ytitle("probability")
xtitle("maternal age + 6 years")
subtitle("THE EFFECT ON PRETERM BIRTH");
graph save FIGA3_ptb.gph,replace;




**********************************************************************************************;
* EXTRA CODING FOR APPENDIX FIGURE A4: NUMBER OF SPECIALIST VISITS BY MATERNAL AGE (1991-1995);
**********************************************************************************************;
* DEFINING THE POPULATION;
for Z in num 1991/1996:
use pnr foeddato koen kom using "D:\Data\WorkData\703201\raw\befZ.dta",clear
\gen yob=year(foeddato)
\gen mob=month(foeddato)
\gen male=koen=="1"
\keep if yob>=1991 & yob<=1995
\save popZ.dta,replace;

use pop1991.dta, clear;
for Z in num 1992/1996: 
\append using popZ.dta;
sort pnr;
keep if pnr!=pnr[_n-1];
rename pnr pnrb;
destring kom,replace;
gen k=1 if kom<=147;
replace k=2 if kom>147 & kom<=189;
for Z in any 155 185: replace k=1 if kom==Z;
replace k=3 if kom==751;
for Z in any 703 705 709 711 713 715 727 737 739 745 601 609 717 721 731 747 749: replace k=5 if kom==Z;
replace k=6 if k==.;
label define muni 1 "CPH" 2 "SUBCPH" 3 "AARH" 5 "SUB AARH" 6 "NOGROUP", modify;
label values k muni;
save pop.dta,replace;
for Z in num 1991/1996: erase popZ.dta;
tab yob;


use pop.dta,clear;
joinby pnrb using "D:\Data\WorkData\703201\raw\fertilitet.dta";
tab yob;
keep k alderm pnrb pnrm pnrf vaegt foeddato yob mob male svlengde;
 tab yob, su(svl);
 drop if svlengde>44;
 drop if svlengde==.;
 gen svl=svlengde;
 gen svlx= svl*7;
 gen concep=foeddato-svlx;
 gen yoc=year(concep);
 gen moc=month(concep);
 drop svl svlx;
  save fertil.dta,replace;

* FINDING THE MOTHERS;
for Z in num 1990/1996:
\use pnrm using fertil.dta,clear
\sort pnrm
\keep if pnrm!=pnrm[_n-1]
\rename pnrm pnr
\joinby pnr using "D:\Data\WorkData\703201\raw\befZ.dta",unmatched(master) 
\tab _m
\drop _m
\gen year=Z
\gen yobm=year(foeddato)
\rename pnr pnrm
\keep pnrm yobm year
\save pnrZ.dta,replace;

use pnr1990.dta, clear;
for Z in num 1991/1996: 
\append using pnrZ.dta;
save popm.dta,replace;
for Z in num 1990/1996: erase pnrZ.dta;



use pnrm using fertil.dta,clear;
joinby pnrm using popm.dta,unmatched(master);
tab _m;
keep if _m==3;
drop if yobm==.;
sort pnrm year;
keep if pnrm!=pnrm[_n+1];
keep pnrm yobm;
save poppnrm.dta,replace;
erase popm.dta;

use poppnrm.dta,clear;
joinby pnrm using "D:\Data\WorkData\703201\raw\fertilitet.dta";
keep pnrm yobm fodtdato  pnrb;
sort pnrm fodtdato;
by pnrm: gen oob=_n;
gen yobb=year(fodtdato);
gen mobb=month(fodtdato);
gen agem =(yobb-yobm);
keep agem yobm oob pnrb pnrm yobb mobb;
tab yobm;
joinby pnrb using fertil.dta;
tab yobm;
keep yobm k oob pnrb pnrm agem yoc moc yobb mobb;
save pop.dta,replace;

use pop.dta,clear;
keep pnrm;
sort pnrm;
keep if pnrm!=pnrm[_n-1];
save poppnrm,replace;
* ADDING INFO FROM THE HEALTH INSURANCE REGISTER;
use "D:\Data\WorkData\703201\raw\syhb1991.dta", clear;
gen year=1991;
for Z in num 1992/1996:
\ append using "D:\Data\WorkData\703201\raw\syhbZ.dta"
\replace year=Z if year==.;
keep year pnr ant*;
rename pnr pnrm;
joinby pnrm using poppnrm.dta;
gen yoc=year;
save sysbh.dta,replace;


use pop.dta,clear;
keep if oob==1;
keep if k==1 | k==2;
keep if yoc==yobb;
tab yoc;
joinby pnrm yoc using sysbh.dta;
keep if yoc<=1995;

gen age=.;
replace age=20 if agem<=20;
replace age=22 if agem>20 & agem<=22;
replace age=24 if agem>22 & agem<=24;
replace age=26 if agem>24 & agem<=26;
replace age=28 if agem>26 & agem<=28;
replace age=30 if agem>28 & agem<=30;
replace age=32 if agem>30 & agem<=32;
replace age=34 if agem>32 & agem<=34;
replace age=36 if agem>34 & agem<=36;
replace age=38 if agem>36 & agem<=38;
replace age=40 if agem>38 ;

label define age  20 "-20" 22 "21-22" 24 "23-24" 26 "25-26" 28 "27-28" 30 "29-30" 32 "31-32" 34 "33-34" 36 "35-36" 38 "37-38" 40 "39+";
label values age age;
drop if antalm>41;
drop if antspec>65;
tab yoc;

**************************************************************************************;
** APPENDIX FIGURE A4.-- THE NUMBER OF GP VISITS AND SPECIALIST VISITS BY MATERNAL AGE;
**************************************************************************************;

graph bar (mean) antalm antspec if oob==1,  legend( label(1 "GP visits") label(2 "Specialist visits (all types)") ) over(age) ytitle("Number of visits") ylabel(0 (2) 12);

save figA4.gph,replace;


**************************************************************************************;
* END OF FIGURES;
**************************************************************************************;


exit, clear;





